module GOLD_surface_forcing
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of GOLD.                                        *
!*                                                                     *
!* GOLD is free software; you can redistribute it and/or modify it and *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* GOLD is distributed in the hope that it will be useful, but WITHOUT *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************
!
!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, November 1998 - May 2002                       *
!*                                                                     *
!*    This program contains the subroutines that calculate the         *
!*  surface wind stresses and fluxes of buoyancy or temperature and    *
!*  fresh water.  These subroutines will be called every time step,    *
!*  even if the wind stresses or buoyancy fluxes are constant in time  *
!*  - in that case these routines return quickly without doing         *
!*  anything.  In addition, any I/O of forcing fields is controlled    *
!*  by surface_forcing_init, located in this file.                     *
!*                                                                     *
!*    set_forcing is a small entry subroutine for the subroutines in   *
!*  this file.  It provides the external access to these subroutines.  *
!*                                                                     *
!*    wind_forcing determines the wind stresses and places them into   *
!*  taux[][] and tauy[][].  Often wind_forcing must be tailored for    *
!*  a particular application - either by specifying file and variable  *
!*  names or by providing appropriate internal expressions for the     *
!*  stresses.                                                          *
!*                                                                     *
!*    buoyancy_forcing determines the surface fluxes of buoyancy,      *
!*  temperature, and fresh water, as is appropriate.  A restoring      *
!*  boundary condition is implemented, but the code for any other      *
!*  boundary condition will usually be modified - either to specify    *
!*  file and variable names and which time level to read, or to set    *
!*  an internal expression for the variables.                          *
!*                                                                     *
!*  Macros written all in capital letters are defined in GOLD_memory.h.*
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v, tauy                                  *
!*    j    x ^ x ^ x   At >:  u, taux                                  *
!*    j    > o > o >   At o:  h, fluxes.                               *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**
!### use GOLD_controlled_forcing, only : apply_ctrl_forcing, register_ctrl_forcing_restarts
!### use GOLD_controlled_forcing, only : controlled_forcing_init, controlled_forcing_end
!### use GOLD_controlled_forcing, only : ctrl_forcing_CS
use GOLD_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use GOLD_cpu_clock, only : CLOCK_MODULE
use GOLD_diag_mediator, only : post_data, query_averaging_enabled
use GOLD_diag_mediator, only : register_diag_field, diag_ptrs, safe_alloc_ptr
use GOLD_domains, only : pass_var, pass_vector, AGRID, To_South, To_West, To_All
use GOLD_error_handler, only : GOLD_error, FATAL, WARNING, GOLD_mesg, is_root_pe
use GOLD_file_parser, only : get_param, log_version, param_file_type
use GOLD_file_parser, only : uppercase
use GOLD_grid, only : ocean_grid_type
use GOLD_initialization, only : Get_GOLD_Input
use GOLD_io, only : file_exists, read_data, slasher, vardesc
use GOLD_restart, only : register_restart_field, restart_init, GOLD_restart_CS
use GOLD_restart, only : restart_init_end, save_restart, restore_state
use GOLD_time_manager, only : time_type, operator(+), operator(/), get_time, set_time
use GOLD_tracer_flow_control, only : call_tracer_set_forcing
use GOLD_tracer_flow_control, only : tracer_flow_control_CS
use GOLD_variables, only : forcing, surface, directories
! use MESO_surface_forcing, only : MESO_wind_forcing, MESO_buoyancy_forcing
! use MESO_surface_forcing, only : MESO_surface_forcing_init, MESO_surface_forcing_CS
use user_surface_forcing, only : USER_wind_forcing, USER_buoyancy_forcing
use user_surface_forcing, only : USER_surface_forcing_init, user_surface_forcing_CS
use user_revise_forcing, only : user_alter_forcing, user_revise_forcing_init
use user_revise_forcing, only : user_revise_forcing_CS
!   Forcing is a structure containing pointers to the forcing fields
! which may be used to drive GOLD.  All fluxes are positive downward.
!   Surface is a structure containing pointers to various fields that
! may be used describe the surface state of GOLD.

implicit none ; private

#include <GOLD_memory.h>

public set_forcing, surface_forcing_init, average_forcing, forcing_save_restart

type, public :: surface_forcing_CS ; private
  logical :: use_temperature ! If true, temperature and salinity are used as
                             ! state variables.
  logical :: restorebuoy     ! If true, use restoring surface buoyancy forcing.
  logical :: adiabatic       ! If true, there are no diapycnal mass fluxes or
                             ! surface buoyancy forcing.
  logical :: variable_winds  ! If true, wind stresses vary with time.
  logical :: variable_buoyforce ! If true, buoyancy forcing varies with time.
  real :: south_lat          ! The southern latitude of the domain.
  real :: len_lat            ! The domain length in latitude.
  real :: Rho0               !   The density used in the Boussinesq
                             ! approximation, in kg m-3.
  real :: G_Earth            !   The gravitational acceleration in m s-2.
  real :: Flux_const         !   The restoring rate at the surface, in m s-1.
  real :: gust_const         !   A constant unresolved background gustiness
                             ! that contributes to ustar, in Pa.
  logical :: read_gust_2d    !   If true, use a 2-dimensional gustiness supplied
                             ! from an input file.
  real, pointer :: gust(:,:) => NULL()    ! A spatially varying unresolved
                             ! background gustiness that contributes to ustar,
                             ! in Pa. gust is used when read_gust_2d is true.
  real, pointer :: T_Restore(:,:) => NULL()  ! The temperature to restore the
                                             ! SST to, in C.
  real, pointer :: S_Restore(:,:) => NULL()  ! The salinity to restore the sea
                                             ! surface salnity to, in PSU.
  real, pointer :: Dens_Restore(:,:) => NULL() ! The density to restore the
                                             ! surface density to, in kg m-3.
  integer :: wind_last_lev_read = -1 ! The last time level read from the
                             ! wind input files.
  integer :: buoy_last_lev_read = -1 ! The last time level read from the
                             ! buoyancy input files.
  real :: gyres_taux_const, gyres_taux_sin_amp, gyres_taux_cos_amp, gyres_taux_n_pis
                             ! if WIND_CONFIG=='gyres' then use
                             ! = A, B, C and n respectively for 
                             ! taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L)
  real :: T_north, T_south   ! Target temperatures at north and south used in
                             ! buoyancy_forcing_linear.
  real :: S_north, S_south   ! Target salinity at north and south used in
                             ! buoyancy_forcing_linear.
  logical :: first_call_set_forcing = .true.
  real :: wind_scale         ! A value by which wind-stresses are scaled, ND.
  character(len=8)   :: wind_stagger
  type(tracer_flow_control_CS), pointer :: tracer_flow_CSp => NULL()
!###  type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL()
  type(GOLD_restart_CS), pointer :: restart_CSp => NULL()
  type(diag_ptrs), pointer :: diag ! A pointer to a structure of shareable
                             ! ocean diagnostic fields and control variables.
  character(len=200) :: inputdir ! The directory where NetCDF input files are.
  character(len=200) :: wind_config ! Indicator for wind forcing type (2gyre, USER, FILE..)
  character(len=200) :: wind_file   ! If wind_config is "file", file to use
  character(len=200) :: buoy_config ! Indicator for buoyancy forcing type
  character(len=200) :: longwavedown_file
  character(len=200) :: longwaveup_file
  character(len=200) :: evaporation_file
  character(len=200) :: sensibleheat_file
  character(len=200) :: shortwaveup_file
  character(len=200) :: shortwavedown_file
  character(len=200) :: snow_file
  character(len=200) :: precip_file
  character(len=200) :: freshdischarge_file
  character(len=200) :: SSTrestore_file
  character(len=200) :: salinityrestore_file
  character(len=80)  :: stress_x_var, stress_y_var

  integer :: id_taux = -1, id_tauy = -1, id_ustar = -1
  integer :: id_PminusE = -1, id_evap = -1, id_precip = -1
  integer :: id_liq_precip = -1, id_froz_precip = -1, id_virt_precip = -1
  integer :: id_liq_runoff = -1, id_froz_runoff = -1
  integer :: id_runoff_hflx = -1, id_calving_hflx = -1
  integer :: id_Net_Heating = -1, id_sw = -1, id_LwLatSens = -1, id_buoy = -1
  integer :: id_LW = -1, id_lat = -1, id_sens = -1
  integer :: id_psurf = -1, id_saltflux = -1, id_TKE_tidal = -1, id_heat_rest=-1

  type(user_revise_forcing_CS), pointer :: urf_CS => NULL()
  type(user_surface_forcing_CS), pointer :: user_forcing_CSp => NULL()
!  type(MESO_surface_forcing_CS), pointer :: MESO_forcing_CSp => NULL()
end type surface_forcing_CS

integer :: id_clock_forcing

contains

subroutine set_forcing(state, fluxes, day_start, day_interval, G, CS)
  type(surface),         intent(inout) :: state
  type(forcing),         intent(inout) :: fluxes
  type(time_type),       intent(in)    :: day_start
  type(time_type),       intent(in)    :: day_interval
  type(ocean_grid_type), intent(inout) :: G
  type(surface_forcing_CS), pointer    :: CS
! This subroutine calls any of the other subroutines in this file
! that are needed to specify the current surface forcing fields.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day_start - Start time of the fluxes.
!  (in)      day_interval - Length of time over which these fluxes
!                           will be applied.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  real :: dt            ! The length of time in seconds over which
                        ! the fluxes will be applied.
  type(time_type) :: day_center  ! The central time of the fluxes.
  integer :: intdt

  call cpu_clock_begin(id_clock_forcing)

  day_center = day_start + day_interval/2
  call get_time(day_interval, intdt)
  dt = real(intdt)

  if (CS%variable_winds .or. CS%first_call_set_forcing) then
    if (trim(CS%wind_config) == "file") then
      call wind_forcing_from_file(state, fluxes, day_center, G, CS)
    elseif (trim(CS%wind_config) == "2gyre") then
      call wind_forcing_2gyre(state, fluxes, day_center, G, CS)
    elseif (trim(CS%wind_config) == "1gyre") then
      call wind_forcing_1gyre(state, fluxes, day_center, G, CS)
    elseif (trim(CS%wind_config) == "gyres") then
      call wind_forcing_gyres(state, fluxes, day_center, G, CS)
    elseif (trim(CS%wind_config) == "zero") then
      call wind_forcing_zero(state, fluxes, day_center, G, CS)
    elseif (trim(CS%wind_config) == "MESO") then
      call GOLD_error(FATAL, "MESO forcing is not available with the ice-shelf"//&
               "version of GOLD_surface_forcing.")
!      call MESO_wind_forcing(state, fluxes, day_center, G, CS%MESO_forcing_CSp)
    elseif (trim(CS%wind_config) == "USER") then
      call USER_wind_forcing(state, fluxes, day_center, G, CS%user_forcing_CSp)
    elseif (CS%variable_winds .and. .not.CS%first_call_set_forcing) then
      call GOLD_error(FATAL, &
       "GOLD_surface_forcing: Variable winds defined with no wind config")
    else
       call GOLD_error(FATAL, &
       "GOLD_surface_forcing:Unrecognized wind config "//trim(CS%wind_config))
    endif
  endif
  if ((CS%variable_buoyforce .or. CS%first_call_set_forcing) .and. &
      (.not.CS%adiabatic)) then
    if (trim(CS%buoy_config) == "file") then
      call buoyancy_forcing_from_files(state, fluxes, day_center, dt, G, CS)
    elseif (trim(CS%buoy_config) == "zero") then
      call buoyancy_forcing_zero(state, fluxes, day_center, dt, G, CS)
    elseif (trim(CS%buoy_config) == "linear") then
      call buoyancy_forcing_linear(state, fluxes, day_center, dt, G, CS)
    elseif (trim(CS%buoy_config) == "MESO") then
      call GOLD_error(FATAL, "MESO forcing is not available with the ice-shelf"//&
               "version of GOLD_surface_forcing.")
!      call MESO_buoyancy_forcing(state, fluxes, day_center, dt, G, CS%MESO_forcing_CSp)
    elseif (trim(CS%buoy_config) == "USER") then
      call USER_buoyancy_forcing(state, fluxes, day_center, dt, G, CS%user_forcing_CSp)
    elseif (trim(CS%buoy_config) == "NONE") then
      call GOLD_mesg("GOLD_surface_forcing: buoyancy forcing has been set to omitted.")
    elseif (CS%variable_buoyforce .and. .not.CS%first_call_set_forcing) then
      call GOLD_error(FATAL, &
       "GOLD_surface_forcing: Variable buoy defined with no buoy config.")
    else
       call GOLD_error(FATAL, &
       "GOLD_surface_forcing: Unrecognized buoy config "//trim(CS%buoy_config))
    endif
  endif
  if (associated(CS%tracer_flow_CSp)) then
    if (.not.associated(fluxes%tr_fluxes)) allocate(fluxes%tr_fluxes)
    call call_tracer_set_forcing(state, fluxes, day_start, day_interval, G, &
                                 CS%tracer_flow_CSp)
  endif

  ! Allow for user-written code to alter the fluxes after all the above
  call user_alter_forcing(state, fluxes, day_center, G, CS%urf_CS)

  CS%first_call_set_forcing = .false.

  call cpu_clock_end(id_clock_forcing)
end subroutine set_forcing


subroutine wind_forcing_zero(state, fluxes, day, G, CS)
  type(surface),            intent(inout) :: state
  type(forcing),            intent(inout) :: fluxes
  type(time_type),          intent(in)    :: day
  type(ocean_grid_type),    intent(inout) :: G
  type(surface_forcing_CS), pointer       :: CS
! This subroutine sets the surface wind stresses.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day - Time of the fluxes.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  real :: PI
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

  if (.not.associated(fluxes%taux)) then
    allocate(fluxes%taux(Isdq:Iedq,jsd:jed)) ; fluxes%taux(:,:) = 0.0
  endif
  if (.not.associated(fluxes%tauy)) then
    allocate(fluxes%tauy(isd:ied,Jsdq:Jedq)) ; fluxes%tauy(:,:) = 0.0
  endif
  if (.not.associated(fluxes%ustar)) then
    allocate(fluxes%ustar(isd:ied,jsd:jed)) ; fluxes%ustar(:,:) = 0.0
  endif
!   Set the steady surface wind stresses, in units of Pa.
  PI = 4.0*atan(1.0)
  do j=js,je ; do I=Isq,Ieq
    fluxes%taux(I,j) = 0.0
  enddo ; enddo
  do J=Jsq,Jeq ; do i=is,ie
    fluxes%tauy(i,J) = 0.0
  enddo ; enddo
  if (CS%read_gust_2d) then
    if (associated(fluxes%ustar)) then ; do j=js,je ; do i=is,ie
      fluxes%ustar(i,j) = sqrt(CS%gust(i,j)/CS%Rho0)
    enddo ; enddo ; endif
  else
    if (associated(fluxes%ustar)) then ; do j=js,je ; do i=is,ie
      fluxes%ustar(i,j) = sqrt(CS%gust_const/CS%Rho0)
    enddo ; enddo ; endif
  endif

end subroutine wind_forcing_zero

subroutine wind_forcing_2gyre(state, fluxes, day, G, CS)
  type(surface),            intent(inout) :: state
  type(forcing),            intent(inout) :: fluxes
  type(time_type),          intent(in)    :: day
  type(ocean_grid_type),    intent(inout) :: G
  type(surface_forcing_CS), pointer       :: CS
! This subroutine sets the surface wind stresses.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day - Time of the fluxes.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  real :: PI
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

  if (.not.associated(fluxes%taux)) then
    allocate(fluxes%taux(Isdq:Iedq,jsd:jed)) ; fluxes%taux(:,:) = 0.0
  endif
  if (.not.associated(fluxes%tauy)) then
    allocate(fluxes%tauy(isd:ied,Jsdq:Jedq)) ; fluxes%tauy(:,:) = 0.0
  endif
  if (.not.associated(fluxes%ustar)) then
    allocate(fluxes%ustar(isd:ied,jsd:jed)) ; fluxes%ustar(:,:) = 0.0
  endif

!   Set the steady surface wind stresses, in units of Pa.
  PI = 4.0*atan(1.0)
  do j=js,je ; do I=Isq,Ieq
    fluxes%taux(I,j) = 0.1*(1.0 - cos(2.0*PI*(G%geolatu(I,j)-CS%South_lat) / &
                                      CS%len_lat))
  enddo ; enddo
  do J=Jsq,Jeq ; do i=is,ie
    fluxes%tauy(i,J) = 0.0
  enddo ; enddo

end subroutine wind_forcing_2gyre

subroutine wind_forcing_1gyre(state, fluxes, day, G, CS)
  type(surface),            intent(inout) :: state
  type(forcing),            intent(inout) :: fluxes
  type(time_type),          intent(in)    :: day
  type(ocean_grid_type),    intent(inout) :: G
  type(surface_forcing_CS), pointer       :: CS
! This subroutine sets the surface wind stresses.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day - Time of the fluxes.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  real :: PI
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

  if (.not.associated(fluxes%taux)) then
    allocate(fluxes%taux(Isdq:Iedq,jsd:jed)) ; fluxes%taux(:,:) = 0.0
  endif
  if (.not.associated(fluxes%tauy)) then
    allocate(fluxes%tauy(isd:ied,Jsdq:Jedq)) ; fluxes%tauy(:,:) = 0.0
  endif
  if (.not.associated(fluxes%ustar)) then
    allocate(fluxes%ustar(isd:ied,jsd:jed)) ; fluxes%ustar(:,:) = 0.0
  endif

!   Set the steady surface wind stresses, in units of Pa.
  PI = 4.0*atan(1.0)
  do j=js,je ; do I=Isq,Ieq
    fluxes%taux(I,j) =-0.2*cos(PI*(G%geolatu(I,j)-CS%South_lat)/CS%len_lat)
  enddo ; enddo
  do J=Jsq,Jeq ; do i=is,ie
    fluxes%tauy(i,J) = 0.0
  enddo ; enddo

end subroutine wind_forcing_1gyre

subroutine wind_forcing_gyres(state, fluxes, day, G, CS)
  type(surface),            intent(inout) :: state
  type(forcing),            intent(inout) :: fluxes
  type(time_type),          intent(in)    :: day
  type(ocean_grid_type),    intent(inout) :: G
  type(surface_forcing_CS), pointer       :: CS
! This subroutine sets the surface wind stresses.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day - Time of the fluxes.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  real :: PI, y
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

  if (.not.associated(fluxes%taux)) then
    allocate(fluxes%taux(Isdq:Iedq,jsd:jed)) ; fluxes%taux(:,:) = 0.0
  endif
  if (.not.associated(fluxes%tauy)) then
    allocate(fluxes%tauy(isd:ied,Jsdq:Jedq)) ; fluxes%tauy(:,:) = 0.0
  endif
  if (.not.associated(fluxes%ustar)) then
    allocate(fluxes%ustar(isd:ied,jsd:jed)) ; fluxes%ustar(:,:) = 0.0
  endif

!   Set the steady surface wind stresses, in units of Pa.
  PI = 4.0*atan(1.0)
  do j=jsd,jed ; do I=Isdq,Iedq
    y = (G%geolatu(I,j)-CS%South_lat)/CS%len_lat
    fluxes%taux(I,j) = CS%gyres_taux_const +                            &
             (   CS%gyres_taux_sin_amp*sin(CS%gyres_taux_n_pis*PI*y)    &
               + CS%gyres_taux_cos_amp*cos(CS%gyres_taux_n_pis*PI*y) )
  enddo ; enddo
  do J=Jsdq,Jedq ; do i=isd,ied
    fluxes%tauy(i,J) = 0.0
  enddo ; enddo

! Set the friction velocities.
  do j=js,je ; do i=is,ie
    fluxes%ustar(i,j) = sqrt(sqrt(0.5*(fluxes%tauy(i,j-1)*fluxes%tauy(i,j-1) + &
      fluxes%tauy(i,j)*fluxes%tauy(i,j) + fluxes%taux(i-1,j)*fluxes%taux(i-1,j) + &
      fluxes%taux(i,j)*fluxes%taux(i,j)))/CS%Rho0 + (CS%gust_const/CS%Rho0))
  enddo ; enddo

end subroutine wind_forcing_gyres

subroutine wind_forcing_from_file(state, fluxes, day, G, CS)
  type(surface),            intent(inout) :: state
  type(forcing),            intent(inout) :: fluxes
  type(time_type),          intent(in)    :: day
  type(ocean_grid_type),    intent(inout) :: G
  type(surface_forcing_CS), pointer       :: CS

! This subroutine sets the surface wind stresses.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day - Time of the fluxes.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  integer :: isd, ied, jsd, jed, Isdq, Iedq, Jsdq, Jedq
  integer :: time_lev = 0          ! With fields from a file, this must
                                   ! be reset, depending on the time.
  character(len=200) :: filename   ! The name of the input file.
  real :: temp_x(SZI_(G),SZJ_(G)) ! Pseudo-zonal and psuedo-meridional
  real :: temp_y(SZI_(G),SZJ_(G)) ! wind stresses at h-points, in Pa.
  integer :: days, seconds

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  Isq = G%Iscq ; Ieq = G%Iecq ; Jsq = G%Jscq ; Jeq = G%Jecq
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
  Isdq = G%Isdq ; Iedq = G%Iedq ; Jsdq = G%Jsdq ; Jedq = G%Jedq

  if (.not.associated(fluxes%taux)) then
    allocate(fluxes%taux(Isdq:Iedq,jsd:jed)) ; fluxes%taux(:,:) = 0.0
  endif
  if (.not.associated(fluxes%tauy)) then
    allocate(fluxes%tauy(isd:ied,Jsdq:Jedq)) ; fluxes%tauy(:,:) = 0.0
  endif
  if (.not.associated(fluxes%ustar)) then
    allocate(fluxes%ustar(isd:ied,jsd:jed)) ; fluxes%ustar(:,:) = 0.0
  endif

  call get_time(day,seconds,days)
  time_lev = days - 365*floor(real(days) / 365.0) +1

  if (time_lev /= CS%wind_last_lev_read) then
    filename = trim(CS%inputdir) // trim(CS%wind_file)
!    if (is_root_pe()) &
!      write(*,'("Wind_forcing Reading time level ",I," last was ",I,".")')&
!           time_lev-1,CS%wind_last_lev_read-1
    select case ( uppercase(CS%wind_stagger(1:1)) )
    case ("A")
      temp_x(:,:) = 0.0 ; temp_y(:,:) = 0.0
      call read_data(filename,CS%stress_x_var,temp_x(:,:), &
                     domain=G%Domain%mpp_domain,timelevel=time_lev)
      call read_data(filename,CS%stress_y_var,temp_y(:,:), &
                     domain=G%Domain%mpp_domain,timelevel=time_lev)

      call pass_vector(temp_x, temp_y, G%Domain, To_All, AGRID)
      do j=js,je ; do I=Isq,Ieq
        fluxes%taux(I,j) = 0.5 * CS%wind_scale * (temp_x(i,j) + temp_x(i+1,j))
      enddo ; enddo
      do J=Jsq,Jeq ; do i=is,ie
        fluxes%tauy(i,J) = 0.5 * CS%wind_scale * (temp_y(i,j) + temp_y(i,j+1))
      enddo ; enddo

      if (CS%read_gust_2d) then
        do j=js,je ; do i=is,ie
          fluxes%ustar(i,j) = sqrt((sqrt(temp_x(i,j)*temp_x(i,j) + &
              temp_y(i,j)*temp_y(i,j)) + CS%gust(i,j)) / CS%Rho0)
        enddo ; enddo
      else
        do j=js,je ; do i=is,ie
          fluxes%ustar(i,j) = sqrt(sqrt(temp_x(i,j)*temp_x(i,j) + &
              temp_y(i,j)*temp_y(i,j))/CS%Rho0 + (CS%gust_const/CS%Rho0))
        enddo ; enddo
      endif
    case ("C")
      call read_data(filename,CS%stress_x_var,fluxes%taux(:,:), &
                     domain=G%Domain%mpp_domain,timelevel=time_lev)
      call read_data(filename,CS%stress_y_var,fluxes%tauy(:,:), &
                     domain=G%Domain%mpp_domain,timelevel=time_lev)
      if (CS%wind_scale /= 1.0) then
        do j=js,je ; do I=Isq,Ieq
          fluxes%taux(I,j) = CS%wind_scale * fluxes%taux(I,j)
        enddo ; enddo
        do J=Jsq,Jeq ; do i=is,ie
          fluxes%tauy(i,J) = CS%wind_scale * fluxes%tauy(i,J)
        enddo ; enddo
      endif

      call pass_vector(fluxes%taux, fluxes%tauy, G%Domain, To_All)
      if (CS%read_gust_2d) then
        do j=js, je ; do i=is, ie
          fluxes%ustar(i,j) = sqrt((sqrt(0.5*((fluxes%tauy(i,j-1)**2 + &
            fluxes%tauy(i,j)**2) + (fluxes%taux(i-1,j)**2 + &
            fluxes%taux(i,j)**2))) + CS%gust(i,j)) / CS%Rho0 )
        enddo ; enddo
      else
        do j=js, je ; do i=is, ie
          fluxes%ustar(i,j) = sqrt(sqrt(0.5*((fluxes%tauy(i,j-1)**2 + &
            fluxes%tauy(i,j)**2) + (fluxes%taux(i-1,j)**2 + &
            fluxes%taux(i,j)**2)))/CS%Rho0 + (CS%gust_const/CS%Rho0))
        enddo ; enddo
      endif
    case default
      call GOLD_error(FATAL, "wind_forcing_from_file: Unrecognized stagger "//&
                      trim(CS%wind_stagger)//" is not 'A' or 'C'.")
    end select
    CS%wind_last_lev_read = time_lev
  endif ! time_lev /= CS%wind_last_lev_read

end subroutine wind_forcing_from_file

subroutine buoyancy_forcing_from_files(state, fluxes, day, dt, G, CS)
  type(surface),         intent(inout) :: state
  type(forcing),         intent(inout) :: fluxes
  type(time_type),       intent(in)    :: day
  real,                  intent(in)    :: dt
  type(ocean_grid_type), intent(inout) :: G
  type(surface_forcing_CS), pointer    :: CS
!    This subroutine specifies the current surface fluxes of buoyancy
!  temperature and fresh water.  It may also be modified to add
!  surface fluxes of user provided tracers.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day - Time of the fluxes.
!  (in)      dt - The amount of time over which the fluxes apply.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.

  real :: rhoXcp ! The mean density times the heat capacity, in J m-3 K-1.
  real :: Irho0  ! The inverse of the Boussinesq density, in m3 kg-1.
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed

  integer :: time_lev              ! With fields from a file, this must
                                   ! be reset, depending on the time.
  integer :: time_lev_monthly      ! With fields from a file, this must
                                   ! be reset, depending on the time.
  integer :: days, seconds
  real, dimension(SZI_(G),SZJ_(G)) :: &
    temp, &       ! A 2-d temporary work array with various units.
    SST_anom, &   ! Instantaneous sea surface temperature anomalies from a
                  ! target (observed) value, in deg C.
    SSS_anom, &   ! Instantaneous sea surface salinity anomalies from a target
                  ! (observed) value, in g kg-1.
    SSS_mean      ! A (mean?) salinity about which to normalize local salinity
                  ! anomalies when calculating restorative precipitation
                  ! anomalies, in g kg-1.

!    Latent heat of vaporization at 15 deg C according to appendix of Gill
  real :: latent_heat_evap=2.4663e6

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed

  if (CS%use_temperature) rhoXcp = CS%Rho0 * fluxes%C_p
  Irho0 = 1.0/CS%Rho0

  if ( CS%use_temperature ) then
    ! Specify the fresh water forcing by setting the following, all in units
    ! of kg m-2 s-1 and positive for water fluxes into the ocean.
    if (.not.associated(fluxes%evap)) then
      allocate(fluxes%evap(isd:ied,jsd:jed))
      fluxes%evap(:,:) = 0.0
    endif
    if (.not.associated(fluxes%liq_precip)) then
      allocate(fluxes%liq_precip(isd:ied,jsd:jed))
      fluxes%liq_precip(:,:) = 0.0
    endif
    if (.not.associated(fluxes%froz_precip)) then
      allocate(fluxes%froz_precip(isd:ied,jsd:jed))
      fluxes%froz_precip(:,:) = 0.0
    endif
    if (.not.associated(fluxes%liq_runoff)) then
      allocate(fluxes%liq_runoff(isd:ied,jsd:jed))
      fluxes%liq_runoff(:,:) = 0.0
    endif
    if (.not.associated(fluxes%froz_runoff)) then
      allocate(fluxes%froz_runoff(isd:ied,jsd:jed))
      fluxes%froz_runoff(:,:) = 0.0
    endif
    if (.not.associated(fluxes%virt_precip)) then
      allocate(fluxes%virt_precip(isd:ied,jsd:jed))
      fluxes%virt_precip(:,:) = 0.0
    endif

    !   Specify the fresh water forcing by setting the following, all in units
    ! of W m-2 and positive for heat fluxes into the ocean.
    if (.not.associated(fluxes%sw)) then
      allocate(fluxes%sw(isd:ied,jsd:jed)) ; fluxes%sw(:,:) = 0.0
    endif
    if (.not.associated(fluxes%lw)) then
      allocate(fluxes%lw(isd:ied,jsd:jed)) ; fluxes%lw(:,:) = 0.0
    endif
    if (.not.associated(fluxes%latent)) then
      allocate(fluxes%latent(isd:ied,jsd:jed)) ; fluxes%latent(:,:) = 0.0
    endif
    if (.not.associated(fluxes%sens)) then
      allocate(fluxes%sens(isd:ied,jsd:jed)) ; fluxes%sens(:,:) = 0.0
    endif
    if (CS%restorebuoy) then
      if (.not.associated(CS%T_Restore)) then
        allocate(CS%T_Restore(isd:ied,jsd:jed))
        CS%T_Restore(:,:) = 0.0
      endif
      if (.not.associated(fluxes%heat_restore)) then
        allocate(fluxes%heat_restore(isd:ied,jsd:jed))
        fluxes%heat_restore(:,:) = 0.0
      endif
      if (.not.associated(CS%S_Restore)) then
        allocate(CS%S_Restore(isd:ied,jsd:jed))
        CS%S_Restore(:,:) = 0.0
      endif
    endif
  else ! CS%use_temperature false.
    if (.not.associated(fluxes%buoy)) then
      allocate(fluxes%buoy(isd:ied,jsd:jed)) ; fluxes%buoy(:,:) = 0.0
    endif
    if (CS%restorebuoy .and. .not.associated(CS%Dens_Restore)) then
      allocate(CS%Dens_Restore(isd:ied,jsd:jed))
      CS%Dens_Restore(:,:) = 0.0
    endif
  endif

  ! Read the file containing the buoyancy forcing.
  call get_time(day,seconds,days)

   time_lev = days - 365*floor(real(days) / 365.0)

  if (time_lev < 31) then ; time_lev_monthly = 0
  else if (time_lev < 59)  then ; time_lev_monthly = 1
  else if (time_lev < 90)  then ; time_lev_monthly = 2
  else if (time_lev < 120) then ; time_lev_monthly = 3
  else if (time_lev < 151) then ; time_lev_monthly = 4
  else if (time_lev < 181) then ; time_lev_monthly = 5
  else if (time_lev < 212) then ; time_lev_monthly = 6
  else if (time_lev < 243) then ; time_lev_monthly = 7
  else if (time_lev < 273) then ; time_lev_monthly = 8
  else if (time_lev < 304) then ; time_lev_monthly = 9
  else if (time_lev < 334) then ; time_lev_monthly = 10
  else ; time_lev_monthly = 11
  endif

  time_lev = time_lev+1
  time_lev_monthly = time_lev_monthly+1

  if (time_lev /= CS%buoy_last_lev_read) then

!   if (is_root_pe()) &
!     write(*,'("buoyancy_forcing : Reading time level ",I3,", last was ",I3,".")')&
!          time_lev,CS%buoy_last_lev_read


    call read_data(trim(CS%inputdir)//trim(CS%longwavedown_file), "lwdn_sfc", &
             fluxes%LW(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev)
    call read_data(trim(CS%inputdir)//trim(CS%longwaveup_file), "lwup_sfc", &
             temp(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev)
    do j=js,je ; do i=is,ie ; fluxes%LW(i,j) = fluxes%LW(i,j) - temp(i,j) ; enddo ; enddo

    call read_data(trim(CS%inputdir)//trim(CS%evaporation_file), "evap", &
             temp(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev)
    do j=js,je ; do i=is,ie
      fluxes%latent(i,j) = -latent_heat_evap*temp(i,j)
      fluxes%evap(i,j) = -temp(i,j)
    enddo ; enddo

    call read_data(trim(CS%inputdir)//trim(CS%sensibleheat_file), "shflx", &
             temp(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev)
    do j=js,je ; do i=is,ie ; fluxes%sens(i,j) = -temp(i,j) ; enddo ; enddo

    call read_data(trim(CS%inputdir)//trim(CS%shortwavedown_file), "swdn_sfc", &
             fluxes%sw(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev)
    call read_data(trim(CS%inputdir)//trim(CS%shortwaveup_file), "swup_sfc", &
             temp(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev)
    do j=js,je ; do i=is,ie
      fluxes%sw(i,j) = fluxes%sw(i,j) - temp(i,j)
    enddo ; enddo

    call read_data(trim(CS%inputdir)//trim(CS%snow_file), "snow", &
             fluxes%froz_precip(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev)
    call read_data(trim(CS%inputdir)//trim(CS%precip_file), "precip", &
             fluxes%liq_precip(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev)
    do j=js,je ; do i=is,ie
      fluxes%liq_precip(i,j) = fluxes%liq_precip(i,j) - fluxes%froz_precip(i,j)
    enddo ; enddo

    call read_data(trim(CS%inputdir)//trim(CS%freshdischarge_file), "disch_w", &
             temp(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev_monthly)
    do j=js,je ; do i=is,ie
      fluxes%liq_runoff(i,j) = temp(i,j)*G%IDXDYh(i,j)
    enddo ; enddo
    call read_data(trim(CS%inputdir)//trim(CS%freshdischarge_file), "disch_s", &
              temp(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev_monthly)
    do j=js,je ; do i=is,ie
      fluxes%froz_runoff(i,j) = temp(i,j)*G%IDXDYh(i,j)
    enddo ; enddo

!     Read the SST and SSS fields for damping.
    if (CS%restorebuoy) then !### .or. associated(CS%ctrl_forcing_CSp)) then
      call read_data(trim(CS%inputdir)//trim(CS%SSTrestore_file), "TEMP", &
               CS%T_Restore(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev_monthly)
      call read_data(trim(CS%inputdir)//trim(CS%salinityrestore_file), "SALT", &
               CS%S_Restore(:,:), domain=G%Domain%mpp_domain, timelevel=time_lev_monthly)
    endif
    CS%buoy_last_lev_read = time_lev

!   Mask out land points.
    do j=js,je ; do i=is,ie
      fluxes%evap(i,j) = fluxes%evap(i,j) * G%hmask(i,j)
      fluxes%liq_precip(i,j)  = fluxes%liq_precip(i,j)  * G%hmask(i,j)
      fluxes%froz_precip(i,j) = fluxes%froz_precip(i,j) * G%hmask(i,j)
      fluxes%liq_runoff(i,j)  = fluxes%liq_runoff(i,j)  * G%hmask(i,j)
      fluxes%froz_runoff(i,j) = fluxes%froz_runoff(i,j) * G%hmask(i,j)
      fluxes%LW(i,j) = fluxes%LW(i,j) * G%hmask(i,j)
      fluxes%latent(i,j) = fluxes%latent(i,j) * G%hmask(i,j)
      fluxes%sens(i,j) = fluxes%sens(i,j) * G%hmask(i,j)
      fluxes%sw(i,j) = fluxes%sw(i,j) * G%hmask(i,j)
    enddo ; enddo
  endif ! time_lev /= CS%buoy_last_lev_read

  if (CS%restorebuoy) then
    if (CS%use_temperature) then
      do j=js,je ; do i=is,ie
        if (G%hmask(i,j) > 0) then
          fluxes%heat_restore(i,j) = G%hmask(i,j) * &
              ((CS%T_Restore(i,j) - state%SST(i,j)) * rhoXcp * CS%Flux_const)
          fluxes%virt_precip(i,j) = - (CS%Rho0*CS%Flux_const) * &
              (CS%S_Restore(i,j) - state%SSS(i,j)) / &
              (0.5*(state%SSS(i,j) + CS%S_Restore(i,j)))
        else
          fluxes%heat_restore(i,j) = 0.0
          fluxes%virt_precip(i,j) = 0.0
        endif
      enddo ; enddo
    else
      do j=js,je ; do i=is,ie
        if (G%hmask(i,j) > 0) then
          fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - state%sfc_density(i,j)) * &
                             (CS%G_Earth*CS%Flux_const/CS%Rho0)
        else
          fluxes%buoy(i,j) = 0.0
        endif
      enddo ; enddo
    endif
  else                                              ! not RESTOREBUOY
    if (.not.CS%use_temperature) then
      call GOLD_error(FATAL, "buoyancy_forcing in GOLD_surface_forcing: "// &
                     "The fluxes need to be defined without RESTOREBUOY.")
    endif
  endif                                             ! end RESTOREBUOY

!### if (associated(CS%ctrl_forcing_CSp)) then
!###   do j=js,je ; do i=is,ie
!###     SST_anom(i,j) = state%SST(i,j) - CS%T_Restore(i,j)
!###     SSS_anom(i,j) = state%SSS(i,j) - CS%S_Restore(i,j)
!###     SSS_mean(i,j) = 0.5*(state%SSS(i,j) + CS%S_Restore(i,j))
!###   enddo ; enddo
!###   call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_restore, &
!###                           fluxes%virt_precip, day, dt, G, CS%ctrl_forcing_CSp)
!### endif

  if (associated(fluxes%p_surf)) then
    do j=js,je ; do i=is,ie
      fluxes%p_surf(i,j) = 0.0
    enddo ; enddo
  endif

end subroutine buoyancy_forcing_from_files

subroutine buoyancy_forcing_zero(state, fluxes, day, dt, G, CS)
  type(surface),         intent(inout) :: state
  type(forcing),         intent(inout) :: fluxes
  type(time_type),       intent(in)    :: day
  real,                  intent(in)    :: dt
  type(ocean_grid_type), intent(in)    :: G
  type(surface_forcing_CS), pointer    :: CS
!    This subroutine specifies the current surface fluxes of buoyancy
!  temperature and fresh water.  It may also be modified to add
!  surface fluxes of user provided tracers.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day - Time of the fluxes.
!  (in)      dt - The amount of time over which the fluxes apply.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  integer :: i, j, is, ie, js, je

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec

  if ( CS%use_temperature ) then
  ! Specify the fresh water forcing by setting the following, all in
  ! units of kg m-2 s-1 and positive for mass fluxes into the ocean.
    if (.not.associated(fluxes%evap)) then
      allocate(fluxes%evap(G%isd:G%ied,G%jsd:G%jed))
      fluxes%evap(:,:) = 0.0
    endif
    if (.not.associated(fluxes%liq_precip)) then
      allocate(fluxes%liq_precip(G%isd:G%ied,G%jsd:G%jed))
      fluxes%liq_precip(:,:) = 0.0
    endif
    if (.not.associated(fluxes%froz_precip)) then
      allocate(fluxes%froz_precip(G%isd:G%ied,G%jsd:G%jed))
      fluxes%froz_precip(:,:) = 0.0
    endif
    if (.not.associated(fluxes%liq_runoff)) then
      allocate(fluxes%liq_runoff(G%isd:G%ied,G%jsd:G%jed))
      fluxes%liq_runoff(:,:) = 0.0
    endif
    if (.not.associated(fluxes%froz_runoff)) then
      allocate(fluxes%froz_runoff(G%isd:G%ied,G%jsd:G%jed))
      fluxes%froz_runoff(:,:) = 0.0
    endif
    if (.not.associated(fluxes%virt_precip)) then
      allocate(fluxes%virt_precip(G%isd:G%ied,G%jsd:G%jed))
      fluxes%virt_precip(:,:) = 0.0
    endif

    !   Specify the heat fluxes by setting the following, all in units
    ! of W m-2 and positive for heat fluxes into the ocean.
    if (.not.associated(fluxes%sw)) then
      allocate(fluxes%sw(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%sw(:,:) = 0.0
    endif
    if (.not.associated(fluxes%lw)) then
      allocate(fluxes%lw(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%lw(:,:) = 0.0
    endif
    if (.not.associated(fluxes%latent)) then
      allocate(fluxes%latent(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%latent(:,:) = 0.0
    endif
    if (.not.associated(fluxes%sens)) then
      allocate(fluxes%sens(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%sens(:,:) = 0.0
    endif
  else
    if (.not.associated(fluxes%buoy)) then
      allocate(fluxes%buoy(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%buoy(:,:) = 0.0
    endif
  endif

  ! This case has no surface buoyancy forcing.

  if (CS%use_temperature) then
    do j=js,je ; do i=is,ie
      fluxes%sw(i,j) = 0.0
      fluxes%lw(i,j) = 0.0
      fluxes%latent(i,j) = 0.0
      fluxes%sens(i,j) = 0.0
      fluxes%liq_precip(i,j) = 0.0
    enddo ; enddo
  else
    do j=js,je ; do i=is,ie
      fluxes%buoy(i,j) = 0.0
    enddo ; enddo
  endif
  if (associated(fluxes%p_surf)) then
    do j=js,je ; do i=is,ie
      fluxes%p_surf(i,j) = 0.0
    enddo ; enddo
  endif

end subroutine buoyancy_forcing_zero

subroutine buoyancy_forcing_linear(state, fluxes, day, dt, G, CS)
  type(surface),         intent(inout) :: state
  type(forcing),         intent(inout) :: fluxes
  type(time_type),       intent(in)    :: day
  real,                  intent(in)    :: dt
  type(ocean_grid_type), intent(in)    :: G
  type(surface_forcing_CS), pointer    :: CS
!    This subroutine specifies the current surface fluxes of buoyancy
!  temperature and fresh water.  It may also be modified to add
!  surface fluxes of user provided tracers.
!
! Arguments: state - A structure containing fields that describe the
!                    surface state of the ocean.
!  (out)     fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      day - Time of the fluxes.
!  (in)      dt - The amount of time over which the fluxes apply.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  real :: y, T_restore, S_restore
  integer :: i, j, is, ie, js, je

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec

  if ( CS%use_temperature ) then
  ! Specify the fresh water forcing by setting the following, all in
  ! units of kg m-2 s-1 and positive for mass fluxes into the ocean.
    if (.not.associated(fluxes%evap)) then
      allocate(fluxes%evap(G%isd:G%ied,G%jsd:G%jed))
      fluxes%evap(:,:) = 0.0
    endif
    if (.not.associated(fluxes%liq_precip)) then
      allocate(fluxes%liq_precip(G%isd:G%ied,G%jsd:G%jed))
      fluxes%liq_precip(:,:) = 0.0
    endif
    if (.not.associated(fluxes%froz_precip)) then
      allocate(fluxes%froz_precip(G%isd:G%ied,G%jsd:G%jed))
      fluxes%froz_precip(:,:) = 0.0
    endif
    if (.not.associated(fluxes%liq_runoff)) then
      allocate(fluxes%liq_runoff(G%isd:G%ied,G%jsd:G%jed))
      fluxes%liq_runoff(:,:) = 0.0
    endif
    if (.not.associated(fluxes%froz_runoff)) then
      allocate(fluxes%froz_runoff(G%isd:G%ied,G%jsd:G%jed))
      fluxes%froz_runoff(:,:) = 0.0
    endif
    if (.not.associated(fluxes%virt_precip)) then
      allocate(fluxes%virt_precip(G%isd:G%ied,G%jsd:G%jed))
      fluxes%virt_precip(:,:) = 0.0
    endif

    !   Specify the heat fluxes by setting the following, all in units
    ! of W m-2 and positive for heat fluxes into the ocean.
    if (.not.associated(fluxes%sw)) then
      allocate(fluxes%sw(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%sw(:,:) = 0.0
    endif
    if (.not.associated(fluxes%lw)) then
      allocate(fluxes%lw(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%lw(:,:) = 0.0
    endif
    if (.not.associated(fluxes%latent)) then
      allocate(fluxes%latent(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%latent(:,:) = 0.0
    endif
    if (.not.associated(fluxes%sens)) then
      allocate(fluxes%sens(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%sens(:,:) = 0.0
    endif
    if (.not.associated(fluxes%heat_restore)) then
      allocate(fluxes%heat_restore(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%heat_restore(:,:) = 0.0
    endif
  else
    if (.not.associated(fluxes%buoy)) then
      allocate(fluxes%buoy(G%isd:G%ied,G%jsd:G%jed)) ; fluxes%buoy(:,:) = 0.0
    endif
  endif

  ! This case has no surface buoyancy forcing.

  if (CS%use_temperature) then
    do j=js,je ; do i=is,ie
      fluxes%sw(i,j) = 0.0
      fluxes%lw(i,j) = 0.0
      fluxes%latent(i,j) = 0.0
      fluxes%sens(i,j) = 0.0
      fluxes%liq_precip(i,j) = 0.0
    enddo ; enddo
  else
    do j=js,je ; do i=is,ie
      fluxes%buoy(i,j) = 0.0
    enddo ; enddo
  endif
  if (associated(fluxes%p_surf)) then
    do j=js,je ; do i=is,ie
      fluxes%p_surf(i,j) = 0.0
    enddo ; enddo
  endif

  if (CS%restorebuoy) then
    if (CS%use_temperature) then
      do j=js,je ; do i=is,ie
        y = (G%geolatu(I,j)-CS%South_lat)/CS%len_lat
        T_restore = CS%T_south + (CS%T_north-CS%T_south)*y
        S_restore = CS%S_south + (CS%S_north-CS%S_south)*y
        if (G%hmask(i,j) > 0) then
          fluxes%heat_restore(i,j) = G%hmask(i,j) * &
              ((T_Restore - state%SST(i,j)) * ((CS%Rho0 * fluxes%C_p) * CS%Flux_const))
          fluxes%virt_precip(i,j) = - (CS%Rho0*CS%Flux_const) * &
              (S_Restore - state%SSS(i,j)) / &
              (0.5*(state%SSS(i,j) + S_Restore))
        else
          fluxes%heat_restore(i,j) = 0.0
          fluxes%virt_precip(i,j) = 0.0
        endif
      enddo ; enddo
    else
      call GOLD_error(FATAL, "buoyancy_forcing_linear in GOLD_surface_forcing: "// &
                     "RESTOREBUOY to linear not written yet.")
     !do j=js,je ; do i=is,ie
     !  if (G%hmask(i,j) > 0) then
     !    fluxes%buoy(i,j) = (CS%Dens_Restore(i,j) - state%sfc_density(i,j)) * &
     !                       (CS%G_Earth*CS%Flux_const/CS%Rho0)
     !  else
     !    fluxes%buoy(i,j) = 0.0
     !  endif
     !enddo ; enddo
    endif
  else                                              ! not RESTOREBUOY
    if (.not.CS%use_temperature) then
      call GOLD_error(FATAL, "buoyancy_forcing_linear in GOLD_surface_forcing: "// &
                     "The fluxes need to be defined without RESTOREBUOY.")
    endif
  endif                                             ! end RESTOREBUOY


end subroutine buoyancy_forcing_linear

subroutine average_forcing(fluxes, dt, G, CS)
  type(forcing),         intent(in) :: fluxes
  real,                  intent(in) :: dt
  type(ocean_grid_type), intent(in) :: G
  type(surface_forcing_CS), pointer    :: CS
!   This subroutine offers forcing fields for time averaging.  These
! fields must first be registered in surface_forcing_init (below).
! This subroutine will typically not be modified, except when new
! forcing fields are added.
!
! Arguments: fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields are unallocated.
!  (in)      dt - The amount of time over which to average.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
  integer :: i,j
  real, dimension(SZI_(G),SZJ_(G)) :: sum

  call cpu_clock_begin(id_clock_forcing)

  if (query_averaging_enabled(CS%diag)) then
    if ((CS%id_taux > 0) .and. ASSOCIATED(fluxes%taux)) &
      call post_data(CS%id_taux, fluxes%taux, CS%diag)
    if ((CS%id_tauy > 0) .and. ASSOCIATED(fluxes%tauy)) &
      call post_data(CS%id_tauy, fluxes%tauy, CS%diag)
    if ((CS%id_ustar > 0) .and. ASSOCIATED(fluxes%ustar)) &
      call post_data(CS%id_ustar, fluxes%ustar, CS%diag)

    if (CS%id_PminusE > 0) then
      sum(:,:) = 0.0
      if (ASSOCIATED(fluxes%liq_precip)) sum(:,:) = sum(:,:)+fluxes%liq_precip(:,:)
      if (ASSOCIATED(fluxes%froz_precip)) sum(:,:) = sum(:,:)+fluxes%froz_precip(:,:)
      if (ASSOCIATED(fluxes%evap)) sum(:,:) = sum(:,:)+fluxes%evap(:,:)
      if (ASSOCIATED(fluxes%liq_runoff)) sum(:,:) = sum(:,:)+fluxes%liq_runoff(:,:)
      if (ASSOCIATED(fluxes%froz_runoff)) sum(:,:) = sum(:,:)+fluxes%froz_runoff(:,:)
      if (ASSOCIATED(fluxes%virt_precip)) sum(:,:) = sum(:,:)+fluxes%virt_precip(:,:)
      call post_data(CS%id_PminusE, sum, CS%diag)
    endif

    if ((CS%id_evap > 0) .and. ASSOCIATED(fluxes%evap)) &
      call post_data(CS%id_evap, fluxes%evap, CS%diag)
    if ((CS%id_precip > 0) .and. ASSOCIATED(fluxes%liq_precip) &
         .and. ASSOCIATED(fluxes%froz_precip)) then
      sum(:,:) = fluxes%liq_precip(:,:) + fluxes%froz_precip(:,:)
      call post_data(CS%id_precip, sum, CS%diag)
    endif

    if ((CS%id_liq_precip > 0) .and. ASSOCIATED(fluxes%liq_precip)) &
      call post_data(CS%id_liq_precip, fluxes%liq_precip, CS%diag)
    if ((CS%id_froz_precip > 0) .and. ASSOCIATED(fluxes%froz_precip)) &
      call post_data(CS%id_froz_precip, fluxes%froz_precip, CS%diag)
    if ((CS%id_virt_precip > 0) .and. ASSOCIATED(fluxes%virt_precip)) &
      call post_data(CS%id_virt_precip, fluxes%virt_precip, CS%diag)
    if ((CS%id_liq_runoff > 0) .and. ASSOCIATED(fluxes%liq_runoff)) &
      call post_data(CS%id_liq_runoff, fluxes%liq_runoff, CS%diag)
    if ((CS%id_froz_runoff > 0) .and. ASSOCIATED(fluxes%froz_runoff)) &
      call post_data(CS%id_froz_runoff, fluxes%froz_runoff, CS%diag)

    if ((CS%id_runoff_hflx > 0) .and. ASSOCIATED(fluxes%runoff_hflx)) &
      call post_data(CS%id_runoff_hflx, fluxes%runoff_hflx, CS%diag)
    if ((CS%id_calving_hflx > 0) .and. ASSOCIATED(fluxes%calving_hflx)) &
      call post_data(CS%id_calving_hflx, fluxes%calving_hflx, CS%diag)

    if (CS%id_Net_Heating > 0) then
      sum(:,:) = 0.0
      if (ASSOCIATED(fluxes%LW)) sum(:,:) = sum(:,:) + fluxes%LW(:,:)
      if (ASSOCIATED(fluxes%latent)) sum(:,:) = sum(:,:) + fluxes%latent(:,:)
      if (ASSOCIATED(fluxes%sens)) sum(:,:) = sum(:,:) + fluxes%sens(:,:)
      if (ASSOCIATED(fluxes%SW)) sum(:,:) = sum(:,:) + fluxes%SW(:,:)
      call post_data(CS%id_Net_Heating, sum, CS%diag)
    endif
    if ((CS%id_LwLatSens > 0) .and. ASSOCIATED(fluxes%lw) .and. &
         ASSOCIATED(fluxes%latent) .and. ASSOCIATED(fluxes%sens)) then
      sum(:,:) = (fluxes%lw(:,:) + fluxes%latent(:,:)) + fluxes%sens(:,:)
      call post_data(CS%id_LwLatSens, sum, CS%diag)
    endif

    if ((CS%id_sw > 0) .and. ASSOCIATED(fluxes%sw)) &
      call post_data(CS%id_sw, fluxes%sw, CS%diag)
    if ((CS%id_LW > 0) .and. ASSOCIATED(fluxes%lw)) &
      call post_data(CS%id_LW, fluxes%lw, CS%diag)
    if ((CS%id_lat > 0) .and. ASSOCIATED(fluxes%latent)) &
      call post_data(CS%id_lat, fluxes%latent, CS%diag)
    if ((CS%id_sens > 0) .and. ASSOCIATED(fluxes%sens)) &
      call post_data(CS%id_sens, fluxes%sens, CS%diag)
    if ((CS%id_heat_rest > 0) .and. ASSOCIATED(fluxes%heat_restore)) &
      call post_data(CS%id_heat_rest, fluxes%heat_restore, CS%diag)

    if ((CS%id_psurf > 0) .and. ASSOCIATED(fluxes%p_surf)) &
      call post_data(CS%id_psurf, fluxes%p_surf, CS%diag)
    if ((CS%id_saltflux > 0) .and. ASSOCIATED(fluxes%salt_flux)) &
      call post_data(CS%id_saltflux, fluxes%salt_flux, CS%diag)
    if ((CS%id_TKE_tidal > 0) .and. ASSOCIATED(fluxes%TKE_tidal)) &
      call post_data(CS%id_TKE_tidal, fluxes%TKE_tidal, CS%diag)

    if ((CS%id_buoy > 0) .and. ASSOCIATED(fluxes%buoy)) &
      call post_data(CS%id_buoy, fluxes%buoy, CS%diag)
  endif

  call cpu_clock_end(id_clock_forcing)
end subroutine average_forcing

subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, &
                                filename_suffix)
  type(surface_forcing_CS),   pointer       :: CS
  type(ocean_grid_type),      intent(inout) :: G
  type(time_type),            intent(in)    :: Time
  character(len=*),           intent(in)    :: directory
  logical,          optional, intent(in)    :: time_stamped
  character(len=*), optional, intent(in)    :: filename_suffix
! Arguments: CS - A pointer to the control structure returned by a previous
!                 call to surface_forcing_init.
!  (in)      G - The ocean's grid structure.
!  (in)      Time - The model time at this call.  This is needed for mpp_write calls.
!  (in, opt) directory - An optional directory into which to write these restart files.
!  (in, opt) time_stamped - If true, the restart file names include
!                           a unique time stamp.  The default is false.
!  (in, opt) filename_suffix - An optional suffix (e.g., a time-stamp) to append
!                              to the restart file names.

  if (.not.associated(CS)) return
  if (.not.associated(CS%restart_CSp)) return

  call save_restart(directory, Time, 1, G, CS%restart_CSp, time_stamped)

end subroutine forcing_save_restart

subroutine surface_forcing_init(Time, G, param_file, diag, CS, tracer_flow_CSp)
  type(time_type),           intent(in) :: Time
  type(ocean_grid_type),     intent(in) :: G
  type(param_file_type),     intent(in) :: param_file
  type(diag_ptrs), target,   intent(in) :: diag
  type(surface_forcing_CS),  pointer    :: CS
  type(tracer_flow_control_CS), pointer :: tracer_flow_CSp
! Arguments: Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure containing pointers to common diagnostic fields.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
!  (in)      tracer_flow_CSp - A pointer to the control structure of the tracer
!                              flow control module.
  type(directories)  :: dirs
  logical            :: new_sim
  type(time_type)    :: Time_frc
  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  character(len=40)  :: mod = "GOLD_surface_forcing" ! This module's name.
  character(len=60)  :: axis_units
  character(len=200) :: filename, gust_file ! The name of the gustiness input file.

  if (associated(CS)) then
    call GOLD_error(WARNING, "surface_forcing_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(CS)

  id_clock_forcing=cpu_clock_id('(Ocean surface forcing)', grain=CLOCK_MODULE)
  call cpu_clock_begin(id_clock_forcing)

  CS%diag => diag
  if (associated(tracer_flow_CSp)) CS%tracer_flow_CSp => tracer_flow_CSp

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mod, version, tagname, "")
  call get_param(param_file, mod, "ENABLE_THERMODYNAMICS", CS%use_temperature, &
                 "If true, Temperature and salinity are used as state \n"//&
                 "variables.", default=.true.)
  call get_param(param_file, mod, "INPUTDIR", CS%inputdir, &
                 "The directory in which all input files are found.", &
                 default=".")
  CS%inputdir = slasher(CS%inputdir)

  call get_param(param_file, mod, "ADIABATIC", CS%adiabatic, &
                 "There are no diapycnal mass fluxes if ADIABATIC is \n"//&
                 "true. This assumes that KD = KDML = 0.0 and that \n"//&
                 "there is no buoyancy forcing, but makes the model \n"//&
                 "faster by eliminating subroutine calls.", default=.false.)
  call get_param(param_file, mod, "VARIABLE_WINDS", CS%variable_winds, &
                 "If true, the winds vary in time after the initialization.", &
                 default=.true.)
  call get_param(param_file, mod, "VARIABLE_BUOYFORCE", CS%variable_buoyforce, &
                 "If true, the buoyancy forcing varies in time after the \n"//&
                 "initialization of the model.", default=.true.)

  call get_param(param_file, mod, "BUOY_CONFIG", CS%buoy_config, &
                 "The character string that indicates how buoyancy forcing \n"//&
                 "is specified. Valid options include (file), (zero), \n"//&
                 "(linear), (USER), and (NONE).", fail_if_missing=.true.)
  if (trim(CS%buoy_config) == "file") then
    call get_param(param_file, mod, "LONGWAVEDOWN_FILE", CS%longwavedown_file, &
                 "The file with the downward longwave heat flux, in \n"//&
                 "variable lwdn_sfc.", fail_if_missing=.true.)
    call get_param(param_file, mod, "LONGWAVEUP_FILE", CS%longwaveup_file, &
                 "The file with the upward longwave heat flux, in \n"//&
                 "variable lwup_sfc.", fail_if_missing=.true.)
    call get_param(param_file, mod, "EVAPORATION_FILE", CS%evaporation_file, &
                 "The file with the evaporative moisture flux, in \n"//&
                 "variable evap.", fail_if_missing=.true.)
    call get_param(param_file, mod, "SENSIBLEHEAT_FILE", CS%sensibleheat_file, &
                 "The file with the sensible heat flux, in \n"//&
                 "variable shflx.", fail_if_missing=.true.)
    call get_param(param_file, mod, "SHORTWAVEUP_FILE", CS%shortwaveup_file, &
                 "The file with the upward shortwave heat flux.", &
                 fail_if_missing=.true.)
    call get_param(param_file, mod, "SHORTWAVEDOWN_FILE", CS%shortwavedown_file, &
                 "The file with the downward shortwave heat flux.", &
                 fail_if_missing=.true.)
    call get_param(param_file, mod, "SNOW_FILE", CS%snow_file, &
                 "The file with the downward frozen precip flux, in \n"//&
                 "variable snow.", fail_if_missing=.true.)
    call get_param(param_file, mod, "PRECIP_FILE", CS%precip_file, &
                 "The file with the downward total precip flux, in \n"//&
                 "variable precip.", fail_if_missing=.true.)
    call get_param(param_file, mod, "FRESHDISCHARGE_FILE", CS%freshdischarge_file, &
                 "The file with the fresh and frozen runoff/calving fluxes, \n"//&
                 "invariables disch_w and disch_s.", fail_if_missing=.true.)
    call get_param(param_file, mod, "SSTRESTORE_FILE", CS%SSTrestore_file, &
                 "The file with the SST toward which to restore in \n"//&
                 "variable TEMP.", fail_if_missing=.true.)
    call get_param(param_file, mod, "SALINITYRESTORE_FILE", CS%salinityrestore_file, &
                 "The file with the surface salinity toward which to \n"//&
                 "restore in variable SALT.", fail_if_missing=.true.)
  endif
  call get_param(param_file, mod, "WIND_CONFIG", CS%wind_config, &
                 "The character string that indicates how wind forcing \n"//&
                 "is specified. Valid options include (file), (2gyre), \n"//&
                 "(1gyre), (gyres), (zero), and (USER).", fail_if_missing=.true.)
  if (trim(CS%wind_config) == "file") then
    call get_param(param_file, mod, "WIND_FILE", CS%wind_file, &
                 "The file in which the wind stresses are found in \n"//&
                 "variables STRESS_X and STRESS_Y.", fail_if_missing=.true.)
    call get_param(param_file, mod, "WINDSTRESS_X_VAR",CS%stress_x_var, &
                 "The name of the x-wind stress variable in WIND_FILE.", &
                 default="STRESS_X")
    call get_param(param_file, mod, "WINDSTRESS_Y_VAR", CS%stress_y_var, &
                 "The name of the y-wind stress variable in WIND_FILE.", &
                 default="STRESS_Y")
    call get_param(param_file, mod, "WINDSTRESS_STAGGER",CS%wind_stagger, &
                 "A character indicating how the wind stress components \n"//&
                 "are staggered in WIND_FILE.  This may be A or C for now.", &
                 default="A")
    call get_param(param_file, mod, "WINDSTRESS_SCALE", CS%wind_scale, &
                 "A value by which the wind stresses in WIND_FILE are rescaled.", &
                 default=1.0, units="nondim")
  endif
  if (trim(CS%wind_config) == "gyres") then
    call get_param(param_file, mod, "TAUX_CONST", CS%gyres_taux_const, &
                 "With the gyres wind_config, the constant offset in the \n"//&
                 "zonal wind stress profile: \n"//&
                 "  A in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
                 units="Pa", default=0.0)
    call get_param(param_file, mod, "TAUX_SIN_AMP",CS%gyres_taux_sin_amp, &
                 "With the gyres wind_config, the sine amplitude in the \n"//&
                 "zonal wind stress profile: \n"//&
                 "  B in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
                 units="Pa", default=0.0)
    call get_param(param_file, mod, "TAUX_COS_AMP",CS%gyres_taux_cos_amp, &
                 "With the gyres wind_config, the cosine amplitude in \n"//&
                 "the zonal wind stress profile: \n"//&
                 "  C in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
                 units="Pa", default=0.0)
    call get_param(param_file, mod, "TAUX_N_PIS",CS%gyres_taux_n_pis, &
                 "With the gyres wind_config, the number of gyres in \n"//&
                 "the zonal wind stress profile: \n"//&
                 "  n in taux = A + B*sin(n*pi*y/L) + C*cos(n*pi*y/L).", &
                 units="nondim", default=0.0)
  endif
  call get_param(param_file, mod, "SOUTHLAT", CS%south_lat, &
                 "The southern latitude of the domain or the equivalent \n"//&
                 "starting value for the y-axis.", units=axis_units, default=0.)
  call get_param(param_file, mod, "LENLAT", CS%len_lat, &
                 "The latitudinal or y-direction length of the domain.", &
                 units=axis_units, fail_if_missing=.true.)
   call get_param(param_file, mod, "RHO_0", CS%Rho0, &
                 "The mean ocean density used with BOUSSINESQ true to \n"//&
                 "calculate accelerations and the mass for conservation \n"//&
                 "properties, or with BOUSSINSEQ false to convert some \n"//&
                 "parameters from vertical units of m to kg m-2.", &
                 units="kg m-3", default=1035.0)
  call get_param(param_file, mod, "RESTOREBUOY", CS%restorebuoy, &
                 "If true, the buoyancy fluxes drive the model back \n"//&
                 "toward some specified surface state with a rate \n"//&
                 "given by FLUXCONST.", default= .false.)
  if (CS%restorebuoy) then
    call get_param(param_file, mod, "FLUXCONST", CS%Flux_const, &
                 "The constant that relates the restoring surface fluxes \n"//&
                 "to the relative surface anomalies (akin to a piston \n"//&
                 "velocity).  Note the non-MKS units.", units="m day-1", &
                 fail_if_missing=.true.)
    ! Convert CS%Flux_const from m day-1 to m s-1.
    CS%Flux_const = CS%Flux_const / 86400.0
    if (trim(CS%buoy_config) == "linear") then
      call get_param(param_file, mod, "SST_NORTH", CS%T_north, &
                 "With buoy_config linear, the sea surface temperature \n"//&
                 "at the northern end of the domain toward which to \n"//&
                 "to restore.", units="deg C", default=0.0)
      call get_param(param_file, mod, "SST_SOUTH", CS%T_south, &
                 "With buoy_config linear, the sea surface temperature \n"//&
                 "at the southern end of the domain toward which to \n"//&
                 "to restore.", units="deg C", default=0.0)
      call get_param(param_file, mod, "SSS_NORTH", CS%S_north, &
                 "With buoy_config linear, the sea surface salinity \n"//&
                 "at the northern end of the domain toward which to \n"//&
                 "to restore.", units="PSU", default=35.0)
      call get_param(param_file, mod, "SSS_SOUTH", CS%S_south, &
                 "With buoy_config linear, the sea surface salinity \n"//&
                 "at the southern end of the domain toward which to \n"//&
                 "to restore.", units="PSU", default=35.0)
    endif
  endif
  call get_param(param_file, mod, "G_EARTH", CS%G_Earth, &
                 "The gravitational acceleration of the Earth.", &
                 units="m s-2", default = 9.80)

  call get_param(param_file, mod, "GUST_CONST", CS%gust_const, &
                 "The background gustiness in the winds.", units="Pa", &
                 default=0.02)
  call get_param(param_file, mod, "READ_GUST_2D", CS%read_gust_2d, &
                 "If true, use a 2-dimensional gustiness supplied from \n"//&
                 "an input file", default=.false.)
  if (CS%read_gust_2d) then
    call get_param(param_file, mod, "GUST_2D_FILE", gust_file, &
                 "The file in which the wind gustiness is found in \n"//&
                 "variable gustiness.", fail_if_missing=.true.)
    call safe_alloc_ptr(CS%gust,G%isd,G%ied,G%jsd,G%jed) ; CS%gust(:,:) = 0.0
    filename = trim(CS%inputdir) // trim(gust_file)
    call read_data(filename,'gustiness',CS%gust,domain=G%domain%mpp_domain, &
                   timelevel=1) ! units should be Pa
  endif
  call get_param(param_file, mod, "AXIS_UNITS", axis_units, default="degrees")

!  All parameter settings are now known.

  if (trim(CS%wind_config) == "USER" .or. trim(CS%buoy_config) == "USER" ) then
    call USER_surface_forcing_init(Time, G, param_file, diag, CS%user_forcing_CSp)
  elseif (trim(CS%wind_config) == "MESO" .or. trim(CS%buoy_config) == "MESO" ) then
    call GOLD_error(FATAL, "MESO forcing is not available with the ice-shelf"//&
               "version of GOLD_surface_forcing.")
!    call MESO_surface_forcing_init(Time, G, param_file, diag, CS%MESO_forcing_CSp)
  endif

  CS%id_taux = register_diag_field('ocean_model', 'taux', G%axesu1, Time, &
        'Zonal Wind Stress', 'Pascal')
  CS%id_tauy = register_diag_field('ocean_model', 'tauy', G%axesv1, Time, &
        'Meridional Wind Stress', 'Pascal')
  CS%id_ustar = register_diag_field('ocean_model', 'ustar', G%axesh1, Time, &
      'Surface friction velocity', 'meter second-1')

  if (CS%use_temperature) then
    CS%id_PminusE = register_diag_field('ocean_model', 'PmE', G%axesh1, Time, &
          'Net fresh water flux (P-E+C+R)', 'kilogram meter-2 second-1')
    CS%id_evap = register_diag_field('ocean_model', 'evap', G%axesh1, Time, &
          'Evaporation at ocean surface (usually negative)', 'kilogram meter-2 second-1')
    CS%id_precip = register_diag_field('ocean_model', 'precip', G%axesh1, Time, &
          'Net (liq.+froz.) precipitation into ocean', 'kilogram meter-2 second-1')
    CS%id_froz_precip = register_diag_field('ocean_model', 'froz_precip', G%axesh1, Time, &
          'Frozen Precipitation into ocean', 'kilogram meter-2 second-1')
    CS%id_liq_precip = register_diag_field('ocean_model', 'liq_precip', G%axesh1, Time, &
          'Liquid Precipitation into ocean', 'kilogram meter-2 second-1')
    CS%id_virt_precip = register_diag_field('ocean_model', 'virt_precip', G%axesh1, Time, &
          'Virtual Precipitation due to salt restoring', 'kilogram meter-2 second-1')
    CS%id_froz_runoff = register_diag_field('ocean_model', 'froz_runoff', G%axesh1, Time, &
          'Frozen runoff (calving) into ocean', 'kilogram meter-2 second-1')
    CS%id_liq_runoff = register_diag_field('ocean_model', 'liq_runoff', G%axesh1, Time, &
          'Liquid runoff (rivers) into ocean', 'kilogram meter-2 second-1')
    CS%id_runoff_hflx = register_diag_field('ocean_model', 'runoff_hflx', G%axesh1, Time, &
          'Heat content of liquid runoff (rivers) into ocean', 'Watt meter-2')
    CS%id_calving_hflx = register_diag_field('ocean_model', 'calving_hflx', G%axesh1, Time, &
          'Heat content of liquid runoff (rivers) into ocean', 'Watt meter-2')

    CS%id_Net_Heating = register_diag_field('ocean_model', 'Net_Heat', G%axesh1, Time, &
          'Net Surface Heating of Ocean', 'Watt meter-2')
    CS%id_sw = register_diag_field('ocean_model', 'SW', G%axesh1, Time, &
        'Shortwave radiation flux into ocean', 'Watt meter-2')
    CS%id_LwLatSens = register_diag_field('ocean_model', 'LwLatSens', G%axesh1, Time, &
          'Combined longwave, latent, and sensible heating', 'Watt meter-2')
    CS%id_lw = register_diag_field('ocean_model', 'LW', G%axesh1, Time, &
        'Longwave radiation flux into ocean', 'Watt meter-2')
    CS%id_lat = register_diag_field('ocean_model', 'latent', G%axesh1, Time, &
        'Latent heat flux into ocean', 'Watt meter-2')
    CS%id_sens = register_diag_field('ocean_model', 'sensible', G%axesh1, Time, &
        'Sensible heat flux into ocean', 'Watt meter-2')
    if (CS%restorebuoy) &
      CS%id_heat_rest = register_diag_field('ocean_model', 'heat_rest', G%axesh1, Time, &
            'Restoring surface heat flux into ocean', 'Watt meter-2')

    CS%id_psurf = register_diag_field('ocean_model', 'p_surf', G%axesh1, Time, &
          'Pressure at ice-ocean or atmosphere-ocean interface', 'Pascal')
    CS%id_saltflux = register_diag_field('ocean_model', 'salt_flux', G%axesh1, Time, &
          'Salt flux into ocean at surface', 'kilogram meter-2 second-1')
    CS%id_TKE_tidal = register_diag_field('ocean_model', 'TKE_tidal', G%axesh1, Time, &
          'Tidal source of BBL mixing', 'Watt meter-2')
  else
    CS%id_buoy = register_diag_field('ocean_model', 'buoy', G%axesh1, Time, &
          'Buoyancy forcing', 'meter2 second-3')
  endif

  ! Set up any restart fields associated with the forcing.
  call restart_init(param_file, CS%restart_CSp, "GOLD_forcing.res")
!###  call register_ctrl_forcing_restarts(G, param_file, CS%ctrl_forcing_CSp, &
!###                                      CS%restart_CSp)
  call restart_init_end(CS%restart_CSp)

  if (associated(CS%restart_CSp)) then
    call Get_GOLD_Input(dirs=dirs)

    new_sim = .false.
    if ((dirs%input_filename(1:1) == 'n') .and. &
        (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true.
    if (.not.new_sim) then
      call restore_state(dirs%input_filename, dirs%restart_input_dir, Time_frc, &
                         G, CS%restart_CSp)
    endif
  endif

!###  call controlled_forcing_init(Time, G, param_file, diag, CS%ctrl_forcing_CSp)

  call user_revise_forcing_init(param_file, CS%urf_CS)

  call cpu_clock_end(id_clock_forcing)
end subroutine surface_forcing_init

subroutine surface_forcing_end(CS, fluxes)
  type(surface_forcing_CS), pointer       :: CS
  type(forcing), optional,  intent(inout) :: fluxes
! Arguments:  CS - A pointer to the control structure returned by a previous
!                  call to surface_forcing_init, it will be deallocated here.
!  (inout)    fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.

  if (present(fluxes)) then
    if (associated(fluxes%taux))        deallocate(fluxes%taux)
    if (associated(fluxes%tauy))        deallocate(fluxes%tauy)
    if (associated(fluxes%ustar))       deallocate(fluxes%ustar)
    if (associated(fluxes%buoy))        deallocate(fluxes%buoy)
    if (associated(fluxes%sw))          deallocate(fluxes%sw)
    if (associated(fluxes%lw))          deallocate(fluxes%lw)
    if (associated(fluxes%latent))      deallocate(fluxes%latent)
    if (associated(fluxes%sens))        deallocate(fluxes%sens)
    if (associated(fluxes%evap))        deallocate(fluxes%evap)
    if (associated(fluxes%liq_precip))  deallocate(fluxes%liq_precip)
    if (associated(fluxes%froz_precip)) deallocate(fluxes%froz_precip)
    if (associated(fluxes%liq_runoff))  deallocate(fluxes%liq_runoff)
    if (associated(fluxes%froz_runoff)) deallocate(fluxes%froz_runoff)
    if (associated(fluxes%virt_precip)) deallocate(fluxes%virt_precip)
    if (associated(fluxes%p_surf))      deallocate(fluxes%p_surf)
    if (associated(fluxes%salt_flux))   deallocate(fluxes%salt_flux)
    if (associated(fluxes%TKE_tidal))   deallocate(fluxes%TKE_tidal)
    if (associated(fluxes%ustar_tidal)) deallocate(fluxes%ustar_tidal)
    ! Deallocate any elements of fluxes%tr_fluxes.
    if (associated(fluxes%tr_fluxes))   deallocate(fluxes%tr_fluxes)
  endif
  
!###  call controlled_forcing_end(CS%ctrl_forcing_CSp)

  if (associated(CS)) deallocate(CS)
  CS => NULL()

end subroutine surface_forcing_end

end module GOLD_surface_forcing
